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ABSTRACT 

We present a study of the galaxy population predicted by hydrodynamical simulations of 
' galaxy clusters. These simulations, which are based on the GADGET-2 Tree+SPH code, include 

, gas cooling, star formation, a detailed treatment of stellar evolution and chemical enrichment, 

as well as SN energy feedback in the form of galactic winds. As such, they can be used to 
extract the spectro-photometric properties of the simulated galaxies, which are identified as 
clumps in the distribution of star particles. Simulations have been carried out for a representa- 
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Q ' live set of 19 cluster-sized halos, having mass A/200 in the range 5 X 10 -1.8 X 10 /i M, 
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All simulations have been performed for two choices of the stellar initial mass function (IMF), 
C/3 ■ namely using a standard Salpeter IMF with power-law index x = 1.35, and a top-heavy IMF 

' with X = 0.95. In general, we find that several of the observational properties of the galaxy 

population in nearby clusters are reproduced fairly well by simulations. A Salpeter IMF is 
, . successful in accounting for the slope and the normalization of the color-magnitude relation 

f\ ' for the bulk of the galaxy population. In contrast, the top-heavy IMF produces too red galax- 

ies, as a consequence of their exceedingly large metallicity. Simulated clusters have a relation 
between mass and optical luminosity which generally agrees with observations, both in nor- 
malization and slope. Also in keeping with observational results, galaxies are generally bluer, 
younger and more star forming in the cluster outskirts. However, we find that our simulated 
clusters have a total number of galaxies which is significantly smaller than the observed one, 
falling short by about a factor 2-3. We have verified that this problem does not have an obvi- 
ous numerical origin, such as lack of mass and force resolution. Finally, the brightest cluster 
galaxies are always predicted to be too massive and too blue, when compared to observations. 
This is due to gas overcooling, which takes place in the core regions of simulated clusters, 
even in the presence of the rather efficient supernova feedback used in our simulations. 
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1 INTRODUCTION 

A key question in the study of the formation and evolution of galax- 
ies concerns the relationship between their observational properties 
and the large-scale cosmological environment. In recent years, a 
flourishing of observational campaign has provided a detailed de- 
scription of the evolution of the galaxy population in clusters. In- 
deed, galaxy clusters play a key role in the characterization of the 
galaxy evolution. Each cluster provides a large sample of galaxies, 
all placed at the same redshift. Furthermore, clusters offer the pos- 



sibility of sampling a variety of environments, from their dense core 
regions, to the outskirts where the properties of the cluster galaxy 
population tends to approach that of the field. 

A diversity of the galaxy population in nearb y clusters, with 
respect to that in the field, was noticed already bv lOemleJ il974) 
and by Dressier ( 1980). Rich clusters were shown to contain a 
higher fraction of bulge-dominated (early type and SO) galax- 
ies, and a corresponding ly lower fraction of star f orming galax- 
ies, than poor systems. iButcher & OemleJ il978h noticed that 
moderately distant clusters {z ~ 0.3) have a galaxy popula- 
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tion which is generally bluer than nearby clusters. This result 
has been subsequently extended by a number of analyses and is 
currently interpreted as high redshift clusters having larger frac- 
tions of star-formin g galaxies than loca l clusters (see Poggiantii 
I2OO4I for a review). lBower et'al] il992h first noticed that early- 
type members of the Coma cluster lie on a tight relation in the 
color-magnitude diagram, a result that has been subsequently con- 
firmed by_ajiumber_of_analyses for extended samples o f nearby 
(e.g., IPrugniel & Simiert 19 96: Andreon 2003; Barrientos et alJ 
l2004lL6pez-Cruz et alj2004 lMcIntosh e t ali2005l and references 
therein) and distant clusters, out to 1 (e.g., 'van Dokkum et al.' 
2000; Blakeslee et al. 2003; Andreon et al. 2004; StrazzuUo et al. 
2006 . and references therein). 

As for the luminosity function (LF hereafter) of cluster galax- 
ies, much work has been done in recent years, with different groups 
reaching different conclusions as for its universality, shape and 
faint-end slop e (e.g.. Driver 2004, for a review). For instance, 
IPopesso et alJ (2006a) have recently analyzed the LF for sample 
of clusters iden tified in the X -ray band within the ROSAT All Sky 
Survey (RASS. |VogeJl992h an d cove red by the Sloan Digital Sky 
Survey (SDSS Stoughton et al. 2002). They found no significant 
cluster to cluster variations of the LF, once calculated within the 
same physical radius (r2oo or '■500)1 with a shape fitted by a double 
Schechter function to account for the upturn at faint magnitudes. 
However, these results are at variance with respect to those from 
other analyses (e.g., Adami et al. 2000). 

Within the widely accepted standard ACDM cosmological 
scenario, galaxies arise from the hierarchical assembly of dark mat- 
ter (DM) halos. The gravitational dynamics of these halos is rela- 
tively simple to describe to high pre cision with modern large super- 
computer simulations (e.g.. lSpringel et al. 2005). However, the ob- 
servational properties of galaxies are determined by the combined 
action of the assembly of DM halos and by the physical processes 
which define the evolution of the cosmic baryons. A complex in- 
terplay between radiative gas cooling, star formation, chemical en- 
richment and release of energy feedback from supernovae (SN) and 
active galactic nuclei (AGN) is expected to determine the properties 
of the stellar population in galaxies. At the same time, the cluster 
environment is expected to play a significant role in altering the 
evolution of galaxies. For instance, ram pressure exerted by the hot 
intra-cluster medium (ICM) can lead to the remov al of a substan - 
tial fraction of the interstellar medium (ISM; <Gunn & Gotllll972h . 
thereby a ffecting galaxy mo r phology, star forma tion and luminos- 
ity (e.g., lAbadi et alj[l999l;lKennev et alJl2004 , and references 
therein). 

In this context, semi-analytical models of galaxy formation 
have been used since several years as a flexible tool to study galaxy 
formation within the c o smological hierarchical framework (e.g., 
Kauffma nn etaP Il9 93'; 'Some rville & Primackl [19991; ICole et alj 
2000; Menci et alJl200Z and references therein). A powerful im- 
plementation of this method is that based on the so-called hy- 
brid approach, which combines N-body simulations, to accurately 
trace the merging history of DM halos, and ser ni-an alytic mod- 
els to describe the physics o f the baryons (e.g., iKauffmann et alj 
ll999h . lSpringel et alJOOOlh applied this method to a DM simula- 
tion of a cluster, with high enough resolution to allow them resolv- 
ing the population of dwarf galaxies. As a result, they found that 
several observational properties (e.g., luminosity function, mass- 
to-light ratio and mo rphological types) are rather well reproduced. 
iDiaferio et alj i200llh applied a semi-analytical model to a DM 
simulation of a large cosmological box, with the aim of perform- 
ing a combined study of kinematics, colors and morphologies for 



both cluster and field galaxies. They concluded that a good agree- 
ment with observations holds for cluster galaxies, while colors and 
star formation rates of field galaxies were shown to evolve more 
rapidly than observed. Casagrande & Diaferio (2006) applied the 
same semi-analytical model to a constrained simulation of the lo- 
cal universe and concluded that significant differences exist be- 
tween the observed and the pr edicted properties of th e large-scale 
distribution of galaxy groups. 'De Lucia et al. ( 2004) incorporated 
in their model also a scheme of metal pro duction to f ollow the 
enrichment of ICM and galaxies (see also ICor3l2006h . Among 
their results, they found that the color-magnitude relation (CMR) is 
mainly driven by metallicity effects, the redd er galaxies on the se- 
quence being on average the more metal rich. Lanzoni et al !'^200j 
applied their semi-analytical model to a set of DM cluster simula- 
tions. They also included a prescription to account for the effect of 
ram-pressure stripping of the ISM as the galaxies move in the hot 
cluster atmosphere, and found it to have only a very little effect on 
the galaxy population. 

A complementary approach to the semi-analytical models is 
represented by using full hydrodynamical simulations, which in- 
clude the processes of gas cooling and star formation. The clear 
advantage of this approach, with respect to semi-analytical mod- 
els, is that galaxy formation can be now described by following in 
detail the evolution of the cosmic baryons while they follow the for- 
mation of the cosmic web. However, the limitation of this approach 
lies in its high computational cost, which prevents it to cover wide 
dynamic ranges and to sample in detail the parameter space de- 
scribing the processes of star formation and feedback. For these 
reasons, describing the process of galaxy formation with a self- 
consistent hydrodynamic approach within the typical cosmological 
environment of ^ 10 Mpc, relevant for galaxy clusters, represents 
a challenging task for sim ulations of the present gen eration. 

In a pioneering paper, iMetzler & Evrardh994h studied the ef- 
fect of including galaxies for the energy feedback and chemical en- 
richment of the ICM. Since these simulations did not have enough 
resolution to identify galaxies, they have been placed b y hand, iden- 
tifying them with the peaks of the initial density field. iFrenk et alJ 
( 1996) used for the first time a radiative simulation of a cluster 
and identified galaxies as concentrations of cooled gas. The aim of 
their study was to compare the dynamics of member galaxies to 
that of DM particles. They concluded that galaxies suffer for a sub- 
stantial dynamical bias, a result which has not been confir med by 
more recent hydrodynamic al simulations (e.g., Faltenbach er et alj 
[2005; "Biviano et al.''200^. Thanks to the ever improving super- 
computing capabilities and efficiency of simulation codes, a num- 
ber of groups have recently completed hydrodynamical simulations 
of galaxy clusters, which have good enough resolution to trace the 
galaxy population with better reliability. Nagai & Rravtsov ( 200^ 
used simulations of eight groups and clusters, performed with an 
adaptive mesh refinement code, including star formation, feedback 
from supernovae and chemical enrichment, to describe the spatial 
distribution of galaxies inside clusters. They found that galaxies are 
more centrally concentrated than DM sub-ha los, with their numbe r 
density profile described by a NFW shape iNavarroetaill 19961) . 
although wi th a smaller concen tration parameter than for the DM 
distribution. lRomeo et alji2005h used SPH simulations of two clus- 
ters, including a similar physics, used spectrophotometric code to 
derive galaxy luminosities in different bands. They analyzed the re- 
sulting color-magnitude relations (CMR, hereafter) and luminosity 
functions, claiming for an overall general agreement with observa- 
tions. 

In this paper, we will present a detailed analysis of the galaxy 
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population for a set of 19 simulated clusters, which span the 
mass range from ~ 5 x 10^^ /i"^Mo to ~ 2 x 10^^ /i"^Mq. 
The simulations have been performed with the Tree-SPH code 
GADGET-2 iS pringel 2005). They include the effect of radiative 
cooli ng, an effective model for st ar formation from a multiphase 
ISM ISpringel & Hemauisi2003'3) . a phenomenological recipe for 
galactic winds, a detailed stellar evolution model, thereby account- 
ing also for li fe-times and metal p roduction from different stellar 
populations I Tomatore et al .12004 Tomatore et al. in preparation). 
These simulations have been carried out also by varying both the 
shape of the initial mass function (IMF, hereafter) and the feedback 
strength. The inclusion of a detailed model of chemical enrichment 
allows us to compute luminosities and colors for galaxies of dif- 
ferent metallicitieSj_byusi^ GALAXEV spectrophotometric 
code <Bruzual & Charloll2003ft . In our analysis we will concentrate 
on the properties of galaxy clusters at 2 = 0. As we shall discuss 
through the paper, several observational trends are reproduced quite 
well by our simulations, although a number of significant discrep- 
ancies are found. For this reason, the aim of our analysis will be 
more that of understanding the directions to improve simulations, 
rather than seeking for a best fitting between model predictions and 
observations. 

The plan of the paper is as follows. In Section 2 we provide 
the general characteristics of the simulated clusters and describe the 
relevant features of the GADGET-2 version used for this analysis. In 
Section 3 we will describe the method of galaxy identification and 
how luminosities in different bands are computed. Section 4 con- 
tains the description of the properties of the simulated galaxy pop- 
ulation and their comparison with observational data. In particular, 
we will discuss the radial distribution of galaxies, the CMR, the 
mass-luminosity ratio, the luminosity function, the star formation 
rate and the color and age gradients. Our main results will be sum- 
marized and discussed in Section 5. We will discuss in an Appendix 
the effects of numerical resolution on the stability of the results of 
our analysis. 



2 THE SIMULATIONS 
2.1 The simulated clusters 

Our set of clusters are identified within nine Lagrangian regions, 
centered around as many main clusters. They were extracted from 
a DM-only simulation with a box size of 479/i~^Mpc of a flat 
ACDM model with Q,m ~ 0.3 for the matter density parameter, 
h — 0.7 for the Hubble constant in units of 100 km s~^Mpc^^, 
(Tg — 0.9 for the r.m.s. fluctuation within a top-hat sphere of 
8 fe~ ^Mpc radius and Q h — 0.04 for the baryon density param- 
eter iYoshida et all200ll) . 

Thanks to the fairly large size chosen for these Lagrangian 
regions, several of them contain other interesting clusters, besides 
the main one. In this way, we end up with 19 clusters with mass 
M200 in the range^ 5 x lO" - 1.8 x 10^^ /I'^Mo, out of which 
4 clusters have M200 > lO^'^ /i'^Mq (see TableQ). Mass resolu- 
tion is increased inside the interestin g regions by using the Zoomed 
Initial Condition (ZIC) technique bv lTormen et alJ 11997). Unper- 
turbed particles positions were placed on a 'glass' ^hite 1996), 
and initial displacements were then assigned according to the Zel- 
dovich approximation (e.g. lShandarin & Zeldovichll989h . Besides 



Table 1. Characteristics of the clusters identified within the simulated re- 
gions aX z = 0. Col. 1: name of the simulated region; Col. 2: name of the 
clusters within each region; Col. 3: value of the total mass, M200, con- 
tained within the radius r200 encompassing an average density 200 times 
larger than the critical cosmic density pc (units of IQi* /i-^M©); Col. 4; 
total number of galaxies, N200 , within r200 ■ having a minimum number of 
32 star particles. 



Region name 


Cluster name 


M200 


^^200 


gl 


gl.a 


12.9 


418 




gl.b 


3.55 


149 




gl.c 


1.39 


51 




gl.d 


0.96 


33 




gl.e 


0.64 


35 


gS 


gS.a 


18.4 


589 




g8.b 


1.02 


42 




gS.c 


0.67 


18 




gS.d 


0.59 


26 




g8.e 


0.54 


21 


g51 


gSl.a 


10.9 


371 


g72 


g72.a 


10.7 


440 




g72.b 


1.55 


60 


g676 


g676.a 


0.89 


23 


g914 


g914a 


0.86 


16 


gl542 


gl542.a 


0.89 


34 


g3344 


g3344.a 


0.97 


29 


g6212 


g6212.a 


0.92 


22 



the low-frequency modes, which were taken from the initial con- 
ditions of the parent simulation, the contribution of the newly sam- 
pled high-frequency modes was also added. The mass resolution 
was progressively degraded in more distant regions, so as to save 
computational resources while still correctly describing the large- 
scale tidal field of the cosmological environment. 

Once initial conditions are created, we split particles in the 
high-resolution region into a DM and a gas component, whose 
mass ratio is set to reproduce the assumed cosmic baryon fraction. 
Instead of placing them on top of each other, in order to avoid spu- 
rious numerical effects, we displace gas and DM particles such that 
the centre of mass of each parent particle is preserved and the fi- 
nal gas and dark matter particle distributions are interleaved by one 
mean particle spacing. In the high-resolution region, the masses of 
the DM and gas particles are set to mDM ~ 1.13 x 10^ h~^MQ 
and mgas = 1.7 x 10* H'^Mq, respectively. The Plummer- 
equivalent softening length for the gravitational force is set to 
epi = 5.0/i~^kpc, kept fixed in physical units from z — 5 to 
z — 0, while being epi — 30.0 /i^^kpc in comoving units at higher 
redshift. 

2.2 The code 

Our simulations are based on an evolution of GADGET-2^ 
(^pririgel et al. 2001; Springel 2005), which includes a de- 
tailed treatment of c hemical enrichment from stellar evolution 
fromatore et aljEo04l ; Tomatore et al., in preparation). GADGET- 
2 is a parallel Tree-l-SPH code with fully adaptive time-stepping, 
which includes an integration scheme which explicitly con- 
serves energy and entropy l,Springel & Hernquist 2Q02.). radiative 
cooling, the effect of a uniform and evolving UV background 
faaardt&Madaiill99 (^. star formation from a multiphase inter- 
stellar medium and a prescription for galactic winds triggered by 



^ We define as the mass contained within a radius encompassing a 
mean density equal to Ape, with pc the critical cosmic density. 
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SN explosions fsee lSpringel & HemguisllFzOOS j for a detailed de- 
scription, SH03 hereafter), and a numerical sc heme to suppress art i- 
ficial viscosity far from the shock regions (see Dolae et al. 
In the original version of the code, energy feedback and global 
metallicity were produced only by SNII under the instantaneous- 
recycling approximation (IRA). 

We have suitably modified the simulation code, so as to cor- 
rectly account for the life-times of different stellar populations, 
to follow metal production from both SNIa and II, while self- 
consistently introducing the depen dence of the cooling function on 
metallicity by using the tables bv ^ Sutherland & Dopital 179931) . A 
detailed description of the implementation of these algorithms will 
be presented in a forthcoming paper (Tornatore et al., in prepara- 
tion), while we provide here a short descriptions of the most rele- 
vant features of the code. 

In order to maintain the general approach of the multiphase 
model by SH03, we assume that stars with masses > 40 Mq ex- 
plode into SNII soon after their formation, thereby promptly re- 
leasing energy and metals. In contrast, we correctly account for 
the lifetime of stars having masses smaller than 40 Mq . The sim- 
ulations that we will discuss here use the lifetimes provided by 
^aeder & Mevnet ( 1989), which have been shown to reproduce the 
abundance pattern in the Milky Way ( ChiaoDini et al. 1997). Within 
the stochastic approach to star formation (SH03), each star particle 
is generated with a mass equal to one third of the mass of its parent 
gas particle. 

Therefore, each star particle is considered as a single stellar 
population (SSP, hereafter), with its own mass, metallicity and red- 
shift of formation. For each SSP we compute both the number of 
stars turning into SNII and la at each time-step and the number 
of stars ending their AGB phase. Then we calculate the amount of 
energy and metals produced by each star particle in a given time in- 
terval, decreasing accordingly the mass of the particle. In this way, 
each star particle is characterized by both its initial mass, assigned 
at the time of its formation, and its final mass, which is updated 
during the evolution. Both SNII and SNIa are assumed to release 
10^^ ergs each, while no energy output is associated to the mass 
loss from AGB stars. The relative number of SNII and SNIa de- 
pends on the choice of the stellar initial mass function (IMF). In 
the following, we will assume for the IMF the power-law shape 
dN/d\ogm oc m^^ . Simulations will be run by assuming the 
Salpeter IMF with s = 1.35 ( Salpet eill95'^ Sa- IMF hereafter) and 
a top-heavy IMF with s = 0.95 ( Ar imoto & Yo shii 1987, TH-IMF 
hereafter). 

The SNIa are associated to bina ry systems whose compo - 
nents are in the 0.8-8 Mq mass range fcreggio & Renzinj||l983h . 
while SNII arise from stars with mass > 8 Mq . In the following, 
we will assume that 10 per cent of stars in the 0.8-8 Mq mass 
range belongs to binary systems, which then produces SNIa. We 
use the analytical fitti ng formulas for stel lar yields of SNIa, SNII 
and PNe provided bv lRecchi etail j200lh. and based on the orig- 
inal nucleosynthesis computations of Nomoto et al 1^19971), usmg 
their W7 model, Wooslev & Weaver (1995) and Renzini & Voli 
( 1981). T he formulation for t he SNIa rate has been calculated as 
in lMatteucci & Recchil i200lh . In the simulations that we present, 
besides H and He, we have followed Fe, O, C, Si, Mg, S. Once 
produced by a star particle, metals are spread over the same num- 
ber of neighbours, 64, used for the SPH computations, also using 
the same kernel. We normalize the IMFs in the mass range O.I-lOO 
Mq. Owing to the uncertainty in modelling yields for very mas- 
sive stars, we take yields to be independent of mass above 40 Mq. 
While any uncertainty in the yields of such massive stars has a neg- 



Table 2. Characteristics of the simulations. Col. 1: simulation name; Col. 
2: IMF slope; Col 3: wind speed, (units of kms~^). 



Name 


IMF slope 




Sa 


1.35 


500 


Sa-NW 


1.35 





TH 


0.95 


500 


TH-SW 


0.95 


1000 



ligible effect for the Salpeter IMF, their accurate description (e.g., 
.Thielemann et am996; Heger & Wooslev. .2002) is required when 
using a top-heavier IMF. 

Our prescription to account for stellar evolution in the sim- 
ulations implies a substantial change of the multiphase "effective 
model" by SH03, which we have suitably modified to account for 
(a) the contribution of the energy reservoir provided by the SN 
which are treated outside the IRA, and (h) the metal-dependence 
of the cool ing function, that we introduce using the tables from 
ISutherland & Dop ita 1 199 3). The resulting density threshold for a 
gas particle to become multiphase, thereby being eligible to un- 
dergo star formation, is fixed to uh =0.1 cm~^ at zero metal- 
licity. According to eq.(23) of SH03, this threshold is inversely 
proportional to the cooling function. Since the latter depends on 
metallicity, we take self-consistently a metallicity-dependent star- 
formation threshold for gas having a non-zero metallicity. 

SH03 also provided a phenomenological description for galac- 
tic winds, which are triggered by SN energy release and whose 
strength is regulated by two parameters. The first one gives the wind 
mass loading according to the relation, Mw ~ rjAI^ , where M, is 
the star formation rate. Following SH03, we assume 77 = 3. The 
second parameter is the wind velocity, v^. For the runs based on 
the Salpeter IMF, we always use — 500 km s^^. For the above 
values of 77 and v^, all the energy from S NII is converted in ki- 
netic energy, as in the original SH03 paper. Ispringel & Hernguisl 
1 2003b) made a study of the star formation history predicted by hy- 
drodinamical simulations which include galactic winds with a simi- 
lar velocity. They concluded that the resulting star fraction at 2 = 
and high-z star formation history are comparable to the observed 
ones. In order to verify the effect of galactic ejecta, the g676 and 
g51 regions are also simulated with the Salpeter IMF, but setting to 
zero the wind velocity (Sa-NW runs). As for the runs based on the 
top-heavy IMF, we will also use vw = 500kms~^ for all clus- 
ters, with the exception of g676 and g5I regions, for which we also 
use Vn, ~ 1000 km (TH-SW runs). In this case, the two wind 
speeds correspond to an energy budget of about 0.4 and 1 .4 times 
the energy provided by SNII. An efficiency larger than unity can 
be justified on the ground of the large uncertainties on the actual 
energy released by SNII explosions. In this perspective, we also 
take the value of the wind velocity as a confidence value. A wind 
velocity vw = 1000 km s^^ is intended to represent an extreme 
feedback case, so that comparing the results with the two values of 
wind velocity allows us to check the effect of a stronger feedback 
on the final properties of the galaxy population. We summarize in 
Table|2|the IMFs and feedback used in our simulations. 
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Figure 1. The number density profile of cluster galaxies. Left panel: the profiles of galaxies of different stellai' mass, averaged over all clusters, for the runs 
with Salpeter IMF (filled symbols). Shown with the solid curve is the average DM profile. All profiles are normalized to the total number density within the 
virial radius. Right panel: The number density of galaxies, brighter than r = —18.5, contained within a given radius, normalized to the total number density 
found within the r200. Filled squares and open circles are by combining all the simulated clusters, for the Salpeter and for the top-heavy IMF, respectively. The 
solid curve is the best-fit King model to the number density profiles of cluster galaxies from the analysis of RASS-SDSS data by Pooesso et al. 1 2()06b), here 
plotted with arbitrary normalization. Errorbars in the simulation profile correspond to Poissonian uncertainties. For reasons of clarity they have been plotted 
only for the Salpeter runs. 



3 ASSIGNING LUMINOSITIES TO GALAXIES 

As a first step, we identify galaxies from the distribution of star par- 
ticles by applying the SKID algorithm'* iStadel 2001). We provide 
a short description of how we applied this algorithm, while a more 
detailed discussion and presentation of tests is provided elsewhere 
(Murante et al. 2006, in preparation; see also|Boraani et al. 2006). 
An overall density field is computed by using the distribution of all 
the particle species, by using a SPH spline-kernel. The star parti- 
cles are then moved along the gradient of the density field in steps 
of r/2, where we assume r ~ 3epi. When a particle begins to 
oscillate inside a sphere of radius r/2, it is stopped. Once all parti- 
cles have been moved, they are grouped using a friends-of-friends 
(FOF) algorithm, with linking length r/2, applied to the new par- 
ticle positions. The binding energy of each group identified in this 
way is then used to remove from the group all star particles which 
are recognized as unbound. All particles in a sphere of radius r, 
centered on the center of mass of the group, are used to compute 
such a gravitational binding energy. Finally, we identify as "bona 
fide" galaxies only those SKlD-groups containing at least 32 star 
particles after the removal of unbound stars. 

Since each star particle is treated as a SSP, with formation 
redshift Zf and metallicity Z,we can assign to it luminosities in 
different bands by resorting to a spectrophotometric code, for the 
appropriate IMF used in the corresponding simulation. 

To this purpose, we have used the outputs of the GALAXEV 
code Bruzual & Chariot (2003) to create a grid of metallicity and 
age values for a SSP of IMq. Luminosities in different bands are 
then assigned to each SSP of this grid. Since two different IMFs 
are used for our simulations, this grid is also computed for both 



^ See |http : / /www-hpcc .astro . Washington . edu/ tools /| 
skid . html 



the lSalDete i* l'l955^ and the top-heavy IMF Note that GALAXEV 
assumes that contributions of different metal species to the total 
metallicity are in solar proportions, while this is not necessarily 
true for the star particles in our simulations. For this reason, we use 
the total metallicity of each star particle (i.e., the sum of the con- 
tributions from the different elements) as input for GALAXEV. We 
have verified that, using instead Iron or Oxygen as proxy for the 
global metallicity, our final results are left essentially unchanged. 
Consistent with the stellar evolution implemented in the simulation 
code, GALAXEV accounts for stellar mass loss. Therefore, we use 
the initial mass of each star particle in the simulations as input to 
GALAXEV to compute the corresponding luminosities. For each 
star particle we interpolate its age and metallicity with the appro- 
priate entries of the grid. Finally, we evaluate the luminosity i*,,^ 
of each star particle, which is treated as a SSP of mass and age 
t, in the i' band by: 

L.At) = ^L^ilMQ) (1) 

In this way, the luminosity in the i/ band of each galaxy is given by 
the sum of the luminosities contributed by each member star parti- 
cle. As a final result, for each galaxy our analysis provides stellar 
mass, mean stellar age, metallicity, star-formation rate (SFR), ab- 
solute magnitudes in the U, B, V, R, I, J, K, bolometric standard 
Johnson bands and in the g, u, r, i, z SLOAN bands. 

We note that GALAXEV accounts for metallicity values in 
the range 0.005-2.6 Zq. While only a negligible number of stars 
have a metallicity below the lower limit of this interval, a sizeable 
number of particles, especially for the top-heavy IMF, are found 
with metallicities exceeding the upper limit. Whenever the particle 
metallicity lies outside the above range, they are set to the value of 
the nearest boundary. 
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4 RESULTS 



4.1 The number-density profile of cluster galaxies 

A well established result from collisionless simulations of galaxy 
clusters is that the radial distribution of ga laxy-sized subhalos 
is less concentrated than that of DM (e .g. lohigna et afj l200Ct 
ISprineel et af]l200ll: lOe Luci a et allbOO^) . and also less concen- 
trated and more extended than the observed radial distrib ution of 
cluster galaxies (e.g. Diemand et al. 2004; Gao et al. 2004). While 
some residual numerical overmerging can still be present at the 
high re solution achieved in D M-only simulations, it has been sug- 
gested jPiemand et al 120041) that overmerging may be physical in 
origin and related to the dissipationless dynamics. The possibility 
to include radiative cooling and star formation in hydrodynamical 
simulations allows one to verify whether the same result holds also 
for the galaxies identified from the star distribution. The general 
conclusion from these analyses is that the radial distribution of sim- 
ulated gala xies is indeed more conce ntrated than that of DM sub- 
halos (e.g., iNagai & Kravtsovll2005h . An intermediate approach, 
based on coupling semi-analytical models of gal axy formation 
with high resolution collisionless simulations (e.g..|Sp rin^e et alJ 
2001'; 'Kravtsov et al.l l2004 iGao et aljEool iDe Lucia et alj|2004 
Lanzoni et al. 2005), confirms that the radial distribution of galax- 
ies is more extended compared to DM. Clearly, the dynamics of 
halo formation in these studies is driven by the collisionless com- 
ponent. Therefore, suitable effective recipes should be included to 
prevent physical overmerging of galaxies within DM halos, so as 
to a chieve agreement wi t h the observed radial galaxy distribution 
(e.g. lSprinsel et all200ll ; lDe Lucia et alJ2004h . 

We show in the left panel of FigureQthe number density pro- 
files of cluster galaxies, after averaging over all the simulated clus- 
ters, compared to the corresponding ave rage DM density profile. 
As discussed bv lNagai & Kravtso\li2005l) . selecting galaxies in hy- 
drodynamical simulations of clusters, which include star formation, 
produces profiles which are generally steeper than those of DM ha- 
los. We confirm here that galaxy profiles become closer to the DM 
profile as more massive galaxies are selected, with objects more 
massive than 10^° /i~^M0 tracing a distribution quite close to the 
DM one. This is just the consequence of the improved capability 
of more massive galaxies to preserve their identity within merg- 
ing halos. Indeed, since galaxies are more concentrated than their 
hosting DM halos, they are able to better survive to disruption and 
merging, thereby providing a better sampling of the underlying DM 
distribution. While this trend is generally consistent with observa- 
tions, a close comparison with data requires assigning luminosities 
to simulated galaxies. For this reason, we also show in the right 
panel of Fig0the number density profiles of galaxies brighter that 
r = —18.5 and compare it to the best fit to the observed profiles 
bv lPopess o et al. ( 2006b), who trace these profiles out to ~ 2r2oo. 
Galaxies of this luminosity have a typical stellar mass of the order 
of 5 X \(f Mq in the simulations with Salpeter IMF. In general, 
the observed and the simulated profiles are quite similar down to 
~ 0.4r2oo. At smaller radii there is a trend for simulated galax- 
ies to have a lower number density than in real clusters. Therefore, 
although tracing galaxies instead of DM halos helps in the com- 
parison with the observed galaxy profiles, still overmerging is also 
the likely reason for the shallower profile as traced by simulated 
galaxies. 





Figure 2. The V-K vs. V color-magnitude relation by combining all the 
galaxies within the virial radii of the simulated clusters, for the Salpeter 
IMF (top panel) and for the top-heavy IMF with normal feedback (bottom 
panel). Straight lin es in each panel show the observed CMR relations by 
iBower et alJ il993i . with the corresponding intrinsic standard deviations. 
Big filled dots mark the BCG of each cluster. Different symbols and colors 
are used for galaxies having different metallicities. Magenta open circles: 
Z > 1.5Zq; blue filled triangles: 1.5 < Z/Zq < 1; red open squares: 
1 < Z/Zq < 0.7; black open triangles: 0.7 < Z/Zq < 0.4; green filled 
squares: Z < 0.4Zq. 




Figure 3. The same as the top panel of Figure|2] but only including in the 
computation of the luminosities the star particles having redshift of forma- 
tion Zf > 1. 



4.2 The color-magnitude relation 

Bright massive ellipticals, which dominate the population of 
cluster galaxies, are observed to form a tight correlation be- 
tween galaxy colors and magnitudes, the so- called red sequence 
or color-magnitud e rela tion (CMR) (e.g.. [ Bower et al] |I992|; 
Prugniel & Simi ei] \l99^ FA ndreon e t all |20()4 iLopez-Cruz et alJ 
2004; Gladders &Yee^'2005t lMcIntoshetalJ l2005^. Attempts to 
compare the observed CMR to that predicted by cosmological mod- 
els of galaxy formation hav e been attempted both using semi- 
analytical approaches (e.g.. iDe Lucia et alj l2004^jLanzoni_etalJ 
2005) and full hydrodynamical simulations lRome^^1l200^ 
As a general results, model predictions reproduce the slope of the 
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CMR reasonably w ell, but with a scatter which is generally larger 
than observed. Rom eo et alj i2005l) simulated one Virgo-sized and 
one Coma-sized cluster. They found that a top-heavy IMF repro- 
duces the normalization of the CMR better than a Salpeter IMF, 
which gives too blue colors as a consequence of the too low metal- 
lic ity. This conclusion i s at variance with respect to that reached 
bv IPe Lucia et alj i2004h who, instead, reproduce the correct CMR 
normalization with a Salpeter IMF. 

Thanks to the larger number of simulated clusters, we can per- 
form the comparison between simulated and observed CMR with 
much improved statistics. The results of this comparison are shown 
in Figure|2| As a term of compariso n, we use the observational de- 
terminations bv lBower et al.|{l992h (see also lTerlevich et all200ll) 
for the V — K \?,.V CMR, which has been determined in the mag- 
nitude range [—18, —23]. Quite apparently, a Salpeter IMF is suc- 
cessful in reproducing the correct amplitude of the CMR at the 
bright end, — 20, while it tends to produce too blue faint 
galaxies. The slope of the CMR appears to be driven by metal- 
licity, the brighter galaxies being redder mainly as a consequence 
of their higher metal conten t, thus in line with the interpretation by 
iKodama & Arimoto At the same time, we note that a top- 

heavy IMF produces too metal-rich galaxies, thereby inducing too 
high a normalization of the CMR. 

The metal content of galaxies is clearly determined by the 
combined action of stellar nucleosynthesis and other processes 
which bring enriched gas far from star forming regions, thus pre- 
venting all metals from being locked back in new ly forming stars. 
Proce sses, such as ram pressure stri pping (e.g.. [Domainko et alJ 
|2005() and galactic winds (e.g., lAguirre et alJl200lh have been 
suggested as the possible mechanisms to enrich the diffuse inter- 
galactic medium. Clearly, the more efficient these mechanisms, the 
lower the expected metallicity of stars and, therefore, the bluer their 
colors. In order to verify whether more efficient galactic winds 
may decrease the metallicity of galaxies in the runs with top- 
heavy IMF, we have re-simulated the g676 and g5 1 clusters using 
Vyj — 1000 km for the galactic winds (TH-SW runs). How- 
ever, while the effect of the stronger feedback is that of decreasing 
the number of galaxies, (see also Sect. l4.4l below'). it leaves their 
metal content, and, therefore, the high CMR normalization, almost 
unchanged. 

Although a Salpeter IMF fares rather well as for the CMR, we 
note that all the BCGs (big filled circles in Fig|2} are much bluer, 
by about 0.5 magnitude, than expected from the red sequence. Such 
a blue excess of the colors of the BGCs, which takes place despite 
their high metallicity, finds its origin in the large star formation 
rate, associated to overcooling, which takes place in the central 
cluster regions. Typical values for the star formation rate of the 
BCG in our simulations are in range 600-1000 A/o/yr for the most 
massive clusters (Af200 10^^ /i'^Mq) and ~ lOOMo/yr for 
the least massive ones (M200 ^ 10^* /i"^M0). Although obser- 
vations indicate the presence of some ongoing star formation in 
some BCGs located at the center of cool core clusters, they are al- 
ways at a much lower level and consistent with a star formation 
rate of ~ lO-lOOAfp/yr for clusters of comparable richness (e.g., 
Ijohnstone et alll987lfBregman et alj2006l:lMcNamara et alj2006l 
and references therein). 

The effect of recent star formation on the CMR is explicitely 
shown in Figure |3] We show here the case in which all star parti- 
cles, formed at redshift z <1 are excluded from the computation of 
the galaxy luminosities. This is equivalent to assume that we com- 
pletely quench star formation since z = 1. Neglecting recent star 
formation has the twofold effect of reducing the scatter in the CMR 



and of making BCG colors significantly redder, although they still 
fall slightly below the observed relation. 

4.3 The mass-luminosity ratio 

A number of observational analyses have established that the 
mass-to-light ratio in clusters generally increases with the clus- 
ter mass, M/L oc IvP with 7 ~ 0.2-0.4, over a fairly 
large dynamic r ange, from poor groups to rich clusters (e.g., 
Adami et allT998 : Girardi et al. 2000, 200 ilBahcall & Comerfor3 
2002; Lin et al. 2003, 2004tlRines et aU2004[lRamella et all2004 
PoDesso et al. 2005). A likely explanation for this trend is the re- 
duced cool ing efficiency within more massive, hotter halos (e.g., 
ISDringel'& HernauisL2003h) . which reduces star formation within 
richer clusters. In fact, an increasing trend of M/L with cluster 
mass is naturally predicted by semi-analytical models of galaxy 
formation (e.g., Kauffmann et al. 1999). 

In Figure|4|we compare the relation between mass and lumi- 
nosity within rso o for our simulated clu sters, and compare it to the 
i- band resul t s bv | PoDesso et alj i2005h and to the if-band results 
bv lLin et al In general, we find that the Af/L from simu- 

lations is rather close to the observed one in the i band, also with a 
comparably small scatter. In the K band, a Salpeter IMF still agrees 
with observations within the statistical uncertainties, while the top- 
heavy IMF produces too red galaxies, thus consistent with the re- 
sults of the CMR, as shown in Fig|2| We fit our mass-luminosity 
relation with a power-law 

1O12L0 ^\10^^MqJ ' ^' 

we find (a,/3), = (0.74,0.92) and (a,/?)^ = (0.76,3.2) 
in the i and K band, respectively, for the runs with Salpeter IMF, 
while (a,/3)i = (0.70,0.91) and {a,(3)K = (0.74,4.7) for the 
top-heavy IMF. Therefore, our simulations agree with the observa- 
tional trend for an increasing mass-to-light ratio with cluster mass, 
independent of the IMF and luminosity band. Using the stronger 
feedback for the top-heavy IMF turns into a sizeable suppression 
of the luminosity, especially for g51. 

The reasonable level of agreement between the observed and 
the simulated A// L may suggest that our simulations produces a re- 
alistic population of galaxies. However, as demonstrated in Figure 
|5| this is not the case. In this figure, we compare the simulated and 
observed number of cluster galaxies, brighter than a given luminos- 
ity limit, both in i and in the K bands. Clearly, simulations under- 
predict such a number, by a factor ~ 2-3. This result is at variance 
with respect to that from semi-analytical models of galaxy forma- 
tion, which instead predict the correct number of cluster members 
(e.g., De Lucia etal. 2004; Lanzoni et al. 2005). However, semi- 
analytical models are generally successful in producing the correct 
LF. They employ a suitable technique to track galaxies, based on 
the assumption that, once a "satellite" galaxy is formed inside a 
DM halo, it preserves its identity and survive to a possible disrup- 
tion of the hosting halo ( Sorineel et al. 2001). Accordingly, the po- 
sition of a galaxy is later assigned to the position of the DM particle 
which was most bound within the DM halo before it was disrupted, 
thereby preventing an excessive merging rate between galaxies. 

On the one hand, it is tempting to explain the lack of galax- 
ies in our simulations as the result of an excessive merging. On the 
other hand, the inclusion of radiative cooling and star formation 
should produce galaxies in our simulations which, in fact, survive 
to the merging of DM halos, and behave as the "satellite" galax- 
ies introduced in the semi-analytical models. Clearly, a reason of 
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Figure 4. The comparison between simulated and observed relation between mass and luminosity in the i band (left panel) and in the K band (right panel). In 
each panel, squares are for the Salpeter IMF, circles for the top-heavy IMF with normal feedback, triangles for the top-heavy IMF with strong feedback and 
asterisks for the Salpeter IMF with no feedback. Filled squares and circles are for the g676 and g5 1 runs, so as to make clear the effect of changing the feedback 
strength. IPopesso et aLi.20Q.5,> for the i band and from .Lin et aL i.2()()41 for the K band, with the dashed lines marking the corresponding observational scatter. 



concern in our simulations is related to the force and mass resolu- 
tions adopted (see Section|2}, which may produce fragile galaxies 
and/or induce spurious numerical overmerging. In the Appendix we 
present a resolution study which is aimed at verifying whether and 
by how much the cluster galaxy population changes when increas- 
ing the resolution. After varying the mass resolution by a factor 45, 
and the corresponding softening parameters by a factor ~ 3.6, we 
find no appreciable variations of the stellar mass function of clus- 
ter galaxies. We also verified that the lack of galaxies is not related 
to numerical heating induced by an non optimal choice of gravita- 
tional softening (e.g., Thomas & Couchman 1992). After running 
a series of simulations, using different choices for epi, we find that 
our softening choice is very close to that maximizing the low end 
of the galaxy stellar mass function. 

As for the effect of feedback, a wind velocity of 500 kms^^ 
is large enough to devoid the gas content of galaxies with mass 
Af^ 10^^ Mq and, therefore, to suppress the number of galax- 
ies above the luminosity limits considered in Fig. |5] Wind veloc- 
ities this high a re generally expected for starburst galaxies (e.g., 
lHeckmaJ200l . while they may be too high for the general galaxy 
population. In order to test this effect, we have performed simu- 
lations of g676 and g5 1 with Salpeter IMF in the extreme case in 
which winds are excluded. In these cases, the numbers of galax- 
ies reported in Fig|5| increase by more than a factor of two, thus 
bringing simulation results into much better agreement with obser- 
vational data. However, the price to pay for this is the increased 
total luminosities, as a result of the larger number of stars formed, 
which introduces a tension between simulations and observations, 
as shown in Fig|3 The need to reconcile at the same time the num- 
ber of galaxies and the total luminosity points toward a scenario in 
which feedback is relatively less effective in small galaxies, while 
being more effective in suppressing star formation in massive rare 
objects. Since massive galaxies are observed to be almost passively 
evolving, this implies that the required feedback mechanism should 
not be directly linked to star formation. In this respect, AGN have 



been suggested to be the natural source fo r this kind of feedback 
('e.g.. lCroton et alJlOOdliower et all2o'o3) . 

4.4 The luminosity function 

The luminosity function (LF hereafter) of cluster galaxies has 
been the subject of numerous studies through the years (e.g. , 
DressleJl97alCollessll989l: lBivian o et alll995l:lGoto et alJ2002t 
D^roprirenLni2003riPop esso et aiJl2006j^, and references 
therein). Despite this, a general consensus on a number of issues 
has still to be reached. Among them, we mention the LF univer- 
sality among clusters and bet ween clusters and field , and the slope 
of the faint end. For instance. IPopesso et alji2006ah have recently 
analysed SDSS data for a set of cluster selected in the X-ray band 
in the RASS. As a result, they found that the LF is universal, once 
calculated within the same physical radiu s, r2oo or rspp. F urther- 
more, the LF can not be fitted by a single ISchechteJ i 19761) func- 
tion, since it displays a marked upturn at faint magnitudes. This 
result is at variance with other analyses. For instance, Adaini et al] 
(2000) performed a deep spectroscopic survey of the Coma cluster 
and found no evidence for an upturn of the LF at faint magnitudes. 

In the following, we will discuss a comparison between the 
LF in our simulated clusters and the observational results by 
Ponesso et al. 1 2006a). To this purpose, we have computed the sim- 
ulated LF, within r2oo, in the r and z b ands, which are two o f the 
four SDSS bands where the analysis bv lPopesso et alj i2006ah has 
been performed. Consist ently with their approach, we have used the 
procedure introduced by Colless I 198?) to compute a composite lu- 
minosity function from the contribution of clusters having different 
richness. Accordingly, the number of galaxies Nj within the j-th 
luminosity bin is defined as 

' ruj ^ No,, 

i 

Here rrij is the number of clusters having galaxies in the j-th lu- 
minosity bin, Nij is the number of galaxies in that luminosity bin 
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Figure 5. The number of galaxies within clusters above a given luminosity limit. Left panel: results in the r band, compared to the observational best-fitting 
result from SDSS data by Pooesso et al. 1 2006b) (the dashed lines ma rk the in trinsic scatter of the observational relation). Right panel: results in the K band, 
compared to the observational best-fitting result from 2MASS data bv lLin et alj 12004) . Symbols for the simulations have the same meaning as in Figurel4l 




Figure 6. The Comparison between the simulated (histograms) and the observed (curves) luminosity functions of cluster galaxies in the Sloan-r (left panel) 
and z (right panel) bands. The smooth curves are the best fit to the SDSS data analysed by Pooesso et al. 1 2006a). In each panel, the solid and the dashed 
histograms are for the Salpeter and for the top-heavy IMF, respectively. Consistent with the observational analysis, the brightest cluster galaxies (BCGs) are 
not included in the computation of the luminosity function. 



contributed by the i-th cluster, A'o,i is the LF normalization for th e 
i-th cluster and A'o = No,i. Followine lPopesso et alj l2006a^ ■ 
we compute A'^o.i as the number of galaxies in the i-th cluster which 
are brighter than r,z — —19. With this definition, each cluster is 
weighted inversely to its richness, in such a way to avoid the rich- 
est clusters to domina te the shape of the LF. Also, consistent with 
IPopesso et alj l2006a h. we do not include the BCGs in the estimate 
of the LF. Owing to the too small number of galaxies found in 
our simulated clusters, we already know in advance that the nor- 
malization of the simulated LF must be lower than the simulated 
one. Therefore, we decide to normalize the simulated LF by hand, 
so that the LF for the Salpeter IMF matches the observed one at 



r = —20 and z = —20. The resulting rescaling factor is then used 
also to re-normalize the LF for the runs with the top-heavy IMF. In 
this way, we preserve the difference in normalization between the 
two series of runs, which is induced by the different choices for the 
IMF. 

The results of this comparison are shown in Figure |6| The 
bright end of the simulated LF is clearly shallower than that of the 
observed one. This is consistent with the picture that overcooling 
takes place within the more massive halos, which hosts the brighter 
galaxies. The simulated LF sh ows a steep ening at the faint end, 
which resembles that found bv lPopesso et a l. (2006a). Looking at 
the combined differential stellar mass function of all the cluster 
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Figure 7. The combined stellar mass function for the galaxies identified 
within r200 of all clusters. The solid and the dashed histograms correspond 
to the Salpeter and to the top-heavy IMF, respectively. 



galaxies (see Figure^, we note an indication for a steepening of 
its slope at the low-mass end, M,^ 3 x 10^'^ Mq. It is this steep- 
ening which causes the corresponding steepening of the luminosity 
functions. Quite interestingly, the faint end of the CMR (see Fig.|2j 
shows a population of small red galaxies, which are in fact associ- 
ated to the excess of faint galaxies shown by the luminosity func- 
tion. It is tempting to make a correspondence between these galax- 
ies eind the faint red galaxies which are claimed by Pooesso et al. 

to contribute to the steepening of their luminosity func- 
tion. However, we consider it as premature to draw strong conclu- 
sions about the slope of the luminosity function in simulations until 
the latter will be demonstrated to roughly produce the correct total 
number of galaxies. 

A comparison between the LPs produced by the Salpeter and 
the top-heavy IMP shows that the latter is shifted towards fainter 
magnitude, especially in the faint end. This effect is also visible in 
the corresponding stellar mass functions of the cluster galaxies (see 
Pig0. While both IMPs produces indistinguishable mass functions 
at the high end, galaxy masses for the top-heavy IMP tend to have 
lower values. This difference is induced by the larger metal con- 
tent associated to the top-heavy IMP, which makes cooling more 
efficient within halos near the resolution limit. Pinally, we show in 
Pigure[8|the effect of increasing the feedback efficiency on the LP. 
In this case, we use the same normalization for the two IMPs, in 
order to directly see the effect of changing the feedback strength. 
Quite interestingly, the effect is that of suppressing the bright end 
of the LP, while leaving the faint end almost unaffected. 

4.5 Radial dependence of the galaxy population 

A number of observations have established that the galaxy pop- 
ulation in clusters is characterized by the presence of color gra- 
dients, with bluer galaxies preferen tially avoiding t o reside in the 
innermost clus t er reg ions (Butcher & Oemle3ll984l . Por instance, 
IPimbblet et alj i200dh found a decreasing trend of the B — R 




Figure 8. The effect of a stronger feedback on the 2-band luminosity func- 
tion. The histograms show the combined luminosity function for the g51 
and g676 clusters in the case of standard feedback {v^ = 500kms~^, 
dashed line) and of strong feedback (v^ = 100 km s~^, solid line). The 
smooth curve is the best-fit to the SDSS data bv lPooesso e7ai]l2006ah . 



color with cluster-centric distance for the galaxies lying on the 
CMR of nearby optically selected cluster s . Similar resul t s have 
als o been f ound by Abraham et alj jl99d) . ICarlberg etall (l9^ 
and lWakeet al. ( 200^ for moderately distant X-ray selected clus- 
ters. Quite consistently, outer cluster regions are populated by a 
larger fraction of blue galaxies (e.g., De Prooris et al. 2004), thus 
confirming that more external galaxies are generally character- 
ized by a relatively younger stellar population. This effect may 
result both as a consequence of the cluster environment, which 
excises star formation in infalling galaxies, and/or due to an ear- 
lier formation ep och of galaxies residing in the cluster center (e.g., 
lEllingson etall20 01). In general, the presence of a gradient in the 
galaxy colors is naturally predicted by semi-analytical models of 
galaxy formation (e.g., Diaferio et al. 2001). 

In Pigure|9|we show the radial variation of the B — V color for 
all galaxies found in our set of simulated clusters. Consistent with 
observational results, the mean galaxy colors become bluer as we 
move towards the outer cluster regions. Quite remarkably, this ef- 
fect extends well beyond the virial radius, thus implying that galax- 
ies feel the cluster environment already at fairly large distances. 
While the trend exists for both a Salpeter and a top-heavy IMP, 
the latter generally predicts much redder colors, consistent with the 
CMR results shown in Pig.|2| Our results for the Salpeter IMP are 
consistent with those reported by Diaferio et al. ( 200 1 ) for the low- 
redshift bin (0.18 <z < 0.3) of the CNOCl cluster sample. 

We note a sudden inversion of the color gradients in the in- 
nermost regions, where galaxies are characterized by much bluer 
colors. These galaxies generally correspond to the cluster BCGs, 
which, as already discussed are much bluer than expected from the 
CMR red sequence (see Pig.|2j. In fact, the galaxies identified in 
this region correspond to the BCG, which, as we have already dis- 
cussed, are characterized by a strong excess of star formation. The 
effect is more pronounced for the top-heavy IMP, which produces 
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Figure 9. The radial dependence of galaxy colors, averaged over all sim- 
ulated clusters. In each panel, solid lines with filled squai'es are for the 
Salpeter IMF, while dashed line with open circles are for the top-heavy 
IMF. The reported results are the average over all the simulated clusters. 
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Figure 10. The fraction of galaxies younger that 8.5 Gyr, as a function of 
cluster-centiic distance, in units of Vyir- Symbols and line types have the 
same meaning as in Figurel9l 



more metals and, therefore, makes overcooling even stronger. This 
highlights once again the presence of too high a cooling rate at 
the center of the simulated clusters, which is not prevented by the 
model of SN feedback adopted in our simulations. 

The presence of color gradients corresponds to the presence of 
age gradients. We show in Figure fTol the radial dependence of the 
fraction of galaxies younger than 8.5 Gyrs. Quite apparently, there 
is a continuous trend for galaxies to be younger in the outer cluster 
regions. The trend extends out to 2rvir, with no evidence for con- 
vergence for a stable mean age in the field. This result further con- 
firms that the presence of a cluster induces environmental effects 
in the galaxy population already at quite large distances. Much like 
for the colors, we note an inversion of the trend in the innermost 
regions, which is due to the excess star formation taking place in 
the central BCGs. The fact that the inversion is more pronounced 
for the top-heavy IMF is in line with its higher enrichment, which 
makes gas cooling more efficient. 

A consistent result also holds for the radial dependence of the 
star formation rate (SFR). In Figure fTTI we show the specific SFR 
(i.e., the SFR per unit stellar mass) as a function of the cluster- 
centric distance. Once we exclude the contribution of the BCG in 
the central bin, we observe a steady increase of the SFR toward ex- 
ternal cluster regions. In general, these results are in line with obser- 
vational evidences for a younger, more star forming galaxy popula- 
tion in the cluster outskirts. For instance, Biviano et al. ( 1997) anal- 
ysed the galaxy population in the ESO Nearby Abell Cluster Survey 
(ENACS) and found that emission-line galaxies te nds to u nderpop- 
ulate the central regions of clusters. Baloeh et al] i 1997b analysed 
data from the CNOCl survey of medium-distant galaxy clusters 
(z ~ 0.2-0.6) and found evidences for a continuous inc rease of the 
SFR out to 2r2oo- In a similar wav. lMoran et alj i2005h analysed a 
large sample of spectroscopic data, covering a 10 Mpc regions 
around a C10024 at z ~ 0.4. Again, they found that galaxies appear 
to be younger at large radii. However, differently from our results. 
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Figure 11. The specific mean star formation rate, averaged over all the sim- 
ulated clusters, as a function of the cluster-centric distance. Symbols and 
line types have the same meaning as in Figure|9| In the central bin we have 
excluded the contribution from the BCGs. 



they detect evidences for an increase of star formation around r^ir, 
possibly triggered by the interaction with the dense environment of 
the ICM. 
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5 CONCLUSIONS 

In this paper we have presented an analysis of the galaxy population 
in cosmological hydrodynamical simulations of galaxy clusters at 
z = Q. The inclusion of a detailed treatment of stellar evolution 
and chemical enrichment iTornatore et al. 2004; Tomatore et al., in 
preparation) in the GADGET-2 code ( Sorinsel 2005) has allowed us 
to derive the properties of galaxies in the optical/near-IR bands. 
In our simulations each star particle is treated as a single stellar 
population (SSP) characterized by a formation time and a metallic- 
ity. Based on this, we apply the GALAXEV spectro-photometric 
code jSmzual & Chariot 2003) to compute luminosities in differ- 
ent bands. Simulations have been carried out for a representative set 
of galaxy clusters, containing 19 objects with mass h'hoo ranging 
from 5 X 10^^ R-^Mq to 1.8 x 10^^ /i^^Mq. All clusters have 
been simulated assuming both a standard Salpeter IMF iSalpeted 
Il955l) and a power-law t op-heavier IMF with exponent x = 0.95 
jArimoto & Yoshiilll987h . The main results of our analysis can be 
summarized as follows. 

1. Both the color-magnitude relation (CMR) and the M/ L ratio 
are in reasonable agreement with observational data for a Salpeter 
IMF. In contrast, using a top-heavy IMF provides too high a metal- 
licity of galaxies, which turns into too red colors. This spoils the 
agreement with the CMR and with the M/L ratio in the K band. 
The CMR is confirmed to be a metallicity sequence, in the sense 
that more enriched galaxies systematically populate the brighter 
redder part of the sequence. 

2. Galaxies are systematically older and redder in central cluster 
regions, thus in keeping with observational results. This trend ex- 
tends at least out to 2rvir, thus showing that the galaxy population 
feels the presence of a cluster well beyond its virial region. Due to 
the overcooling occurring in the central cluster regions, BCGs are 
always much bluer and more massive than observed, and charac- 
terized by too high a recent star formation. Indeed, neglecting the 
contribution of stars formed at 2 < 1 produces significantly redder 
BCGs, thus in better agreement with observational data. 

3. The number density profile of galaxies is confirmed to steepen 
with the galaxy stellar mass, approaching the DM profile for the 
galaxies with M, > 10^° h~^Mp). However, when compared to 
the observations fe.g.. |PoDesso et all2006b.) . the simulation profiles 
are flatter that the observed ones inside 0.4r2oo- 

4. Simulated clusters have about three times fewer galaxies above 
a given luminosity limit than real clusters. We have verified that this 
disagreement is directly related neither to lack of mass and force 
resolution, nor to numerical gas heating due to a non optimal choice 
of the softening parameter. This leaves as further possibilities more 
subtle numerical effects or a better suited implementation of energy 
feedback. 

5. The luminosity function (LF) is shallower than the observed 
one in the bright end, thus confirming that feedback is not strong 
enough to suppress cooling in the most massive halos. In the faint 
end, the LF steepens and indicates the presence of an excess of 
small red galaxies. A lthough this result res embles that found in 
observational data bv lPopesso et alj i2006ah . we warn against its 
overinterpretation, in the light of the deficit of the overall number 
of galaxies found in the simulated clusters. 

The results of our analysis support the capability of hydro- 
dynamical simulations of galaxy clusters to reproduce the general 
trends characterizing their galaxy population. Therefore, simula- 
tions in which the properties of the stellar population are self- 
consistently predicted from the gas dynamics, can provide an alter- 



native and complementary approach to semi-analytical methods. 
However, at present, our simulations have two main limitations in 
accounting for the observed properties of the cluster galaxy popu- 
lation. 

First of all, the deficit of galaxies produced in our simula- 
tions suggests that they are not able to produce galaxies which are 
resistant enough to survive the tidal field of the cluster environ- 
ment. Indeed, the shallow galaxy number density profile in cen- 
tral cluster regions shows that the problem is more apparent where 
effects of galaxy disruption are expected to be stronger. This re- 
sult on the lack of galaxies is at variance with respect to the pre- 
dictions of semi-analytical models (SAM) of galaxy formation, 
which produce rough l y the correct nu mber of cluster galaxies (e.g., 
iDe Lucia et aljE004l : iLanzoni et alJEoOSV Comparisons between 
the galaxy populations predicted by hydrodynamical simulations 
and by SAM have shown a reasonable level of agreement. However, 
these comparisons have been always performed by excluding the 
effect of energy feedback (e.g., iHellv et alj|2003l : ICattaneo et alj 
120061) or explic it conversion of coo led gas particles into collision- 
less stars (e.g., lYoshida et all200^ . 

Furthermore, we always find that the brightest cluster galax- 
ies (BCGs) are much bluer and star forming that observed. While 
the adopted feedback scheme, based on galactic winds, is efficient 
in regulating star formation for the bulk of the galaxy population, 
it is not able to quench low-redshift star formation in the central 
cluster regions, to a level consistent with observations. Clearly, the 
required feedback mechanism should be such to leave the bulk of 
the galaxy population unaffected while acting only on the very 
high end of the galaxy mass distribution. Feedback from central 
AGN represents the natural solu tion and provides i n principle a 
large enough energy budget (e.g., iRaffertv e t al. 2006). Its detailed 
implementation in cosmological hydrodynamical simulations re- 
quire understanding in detail the mechanisms for the thermaliza- 
tion of this energy in the diffuse medium (e.g., IZanni et all2005t 
Siiacki & Sorinael 2006). 

An ambitious goal for hydrodynamical simulations of the next 
generation will be that of describing in detail the complex interplay 
between the history of star formation and the thermodynamical and 
chemical evolution of the diffuse cosmic baryons. For this to be 
reached, it is mandatory that simulations are able to produce a real- 
istic population of galaxies, both in terms of color and of luminosity 
and mass distribution. Therefore, simulation codes will be required 
on the one hand to have numerical effects under exquisite control 
and, on the other hand, to include physically motivated schemes of 
energy and metal feedback. 
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APPENDIX. TESTING NUMERICAL EFFECTS ON THE 
GALAXY STELLAR MASS FUNCTION 

In this Appendix we discuss the stability, against possible numer- 
ical effects, of the stellar mass function of the galaxies identified 
inside simulated clusters. In particular, we will focus the analysis 
on fa) the effect of changing the softening of the gravitational force, 
to control the possible presence of spurious numerical heating (e.g., 
iThomas & Couchman 1992); (b) the effect of mass and force reso- 
lution. As for the softening choice, it is known that using too small 
values may induce spurious heating of the gas particles by two- 
body collisions, thereby inhibiting gas cooling inside small halos. 
On the other hand, increasing it to too large a value also reduces 
the number of galaxi es as a consequence of the lower number of 
resolved small halos <Borgani et alll2006h . As for the resolution, 
increasing it has the effect to better resolve the low end of the mass 
function and, in general, is expected to produce a more reliable 
galaxy population. 

The tests described in this Appendix are based on three clus- 
ter sized halos, which have been resimulated by varying mass 
and force resolution. These clusters have virial masses in the 
range (1.6-2.9) X lO^^/i-^M© and are described in detail by 
■Borgani et aL 12006) . They have been simulated for the same cos- 
mological model of the clusters described in this paper, but with 
a lower normalization of the power spectrum, erg = 0.8. At the 
lowest resolution, the mass of the gas particle is rrigas — 6.9 x 
10* /i~^M0 and the Plummer-equivalent softening for gravita- 
tional forces is set to e = 7.5^~^kpc at z = 0. We point out 
that the simulations analysed in this paper and described in Section 
I2.1l have a mass resolution which is better than the above one by 
about a factor of four. Runs at increasingly higher resolution have 
been performed by decreasing the particle masses by a factor 3, 
10 and 45, with the softening correspondingly decreased accord- 
ing to the m}^^ scaling. Therefore, at the highest resolution, it is 
rrigas — 1.5 X lO^/r~^M0 and e = 2.1/i~^kpc. This mass res- 
olution is ~ 11 times higher than used for the set of simulations 
analysed in this paper. For the most massive of these three clusters, 
we have also repeated the simulation at 10 times the basic mass res- 
olution with four different choices of the gravitational softening. In 
particular, we have decreased it by a factor two, with respect to 
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Figure 12. The cumulative stellar mass function of the galaxies identi- 
fied within the virial radii of the most massive among the three simulated 
clusters described in the Appendix (see alsopSorgani et al. 2006). All the 
simulations have been done at fixed mass resolution, which correspond 
to an increase by a factor 10 with respect to the basic resolution (i.e., 
see text). The four curves correspond to the 
different choices for the gravitational softening. The labels indicate the fac- 
tor by which the softening has been changed, with respect to the standard 
choice of 3.5/i~^kpc. 



the standard choice, and increased it by a factor two and four. As 
such, this set of simulations allows us to verify the stability of the 
galaxy stellar mass function against changing mass resolution and 
the choice for the gravitational softening. 

The simulations have been performed with the origi- 
nal prescription f or sta r forma tion and feedback presented 
by [Sorinsel & He mauis3 J2003ah . with a wind speed v^u — 
480 kms^^, therefore comparable to that assumed for the standard 
feedback in this paper. However, those runs did not include the pre- 
scription for stellar evolution and chemical enrichment, which we 
used here. Since each gas particle is allowed to spawn two star par- 
ticles, the latter have a m ass which is half of t hat of the parent gas 
particle. As discussed by Borsani et al. (2006), this series of runs 
produces an amount of stars within the virial radius of the clusters, 
which is almost independent of the resolution, thereby preventing 
the runaway of cooling with increasing resolution. 

We show in Figure [121 the effect of varying the softening on 
the cumulative stellar mass function of the galaxies identified in- 
side the virial radius. As expected, decreasing the softening to half 
the standard value has the effect of suppressing the low end of the 
mass function, Afi$ 2 x IO^^/i^^A/q, as a consequence of spurious 
numerical heating of gas. On the other hand, increasing the soften- 
ing by a factor four also induces a suppression of low-mass galax- 
ies, as a consequence of the lack of resolution. At larger masses, 
using a too small softening has the effect of increasing the mass 
function, although the rather small number of galaxies in the high 
mass end prevents from detecting systematic trends. These results 
demonstrate that our lack of galaxies can not be explained by a 
non-optimal choice of the gravitational force softening. 

As for the effect of resolution, we plot in Figure com- 
bined cumulative stellar mass function for all the galaxies identi- 
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Figure 13. The combined cumulative stellar mass function of the galaxies 
identified within the virial radii of the three clusters. The four curves corre- 
spond to the different resolutions at which the clusters have been simulated. 
Continuous, long-dashed, short-dashed and dotted curves are for the simu- 
lations at progressively increasing resolution. The labels indicate the factor 
by which mass resolution is increased, with respect to the lowest resolution 
run (Ix). 

fied within the virial radii of the three clusters, simulated at four 
different resolutions. The first apparent effect of increasing res- 
olution is that of steepening the mass function in the low mass 
end. In the mass range where galaxies are identified with at least 
64 star particles at the different resolutions, which correspond to 
Mt ~ 2.2 X 10^°h~'^MQ for the lowest resolution run, the mass 
functions have a weaker dependence on resolution, with a decreas- 
ing trend of the high end of the mass function. This steepening of 
the high end of the mass function at increasing resolution is the con- 
sequence of the reduction of overmerging, which makes small halos 
surviving more efficiently and, therefore, prevents their disruption 
and accretion inside massive halos. This result demonstrates that, 
at least at the highest resolution reached in this test, also resolution 
is not the reason for the too low number of galaxies found in the 
simulated clusters, when compared to observations (see discussion 
in Sect.l43l. 



